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Bilayer honeycomb lattice, with inter-layer tunneling energy, has a parabolic dispersion relation, and the 
inter-layer hopping can cause the charge imbalance between two sublattices. Here, we investigate the 
metal-insulator and magnetic phase transitions on the strongly correlated bilayer honeycomb lattice by 
cellular dynamical mean-field theory combined with continuous time quantum Monte Carlo method. The 
procedures of magnetic spontaneous symmetry breaking on dimer and non-dimer sites are different, 
causing a novel phase transition between normal anti-ferromagnet and layer anti-ferromagnet. The whole 
phase diagrams about the magnetism, temperature, interaction and inter-layer hopping are obtained. 
Finally, we propose an experimental protocol to observe these phenomena in future optical lattice 
experiments. 

The bilayer honeycomb lattice (BHL) has attracted enormous interest in both experimental and theoretical 
research. Lots of novel phenomena have been found in BHL, for instance, the quantum Hall effect, quantum 
spin Hall effect, and chiral superconductivity^"^ \ However, the real charge and magnetic order induced by 
Coulomb interaction are still challenges in the strongly correlated BHL^^"^^. A quadratic dispersion relation, 
signed by two touching bands in the corners of Brillouin zone, can be found in BHL, which is driven by the inter- 
layer hopping. In addition, the spontaneous symmetry can be broken by the dimers when the inter-layer hopping 
changes. Some amazing phases may emerge, such as layer anti-ferromagnetic phase and paramagnetic insulator 
phase. Previous work mainly focus on the electronic properties of BHL^^"^^, and the phase diagram for the 
magnetic phase transition induced by the interaction and dimers are still absent. Moreover, the progress of 
optical lattice provides us a useful tool to set a controllable and clearness experimental platform to simulate 
the strongly correlated BHL, in which the interaction between trapped fermionic cold atoms can be tuned by the 
Feshbach resonance^^"^\ 

The dynamical mean-field theory (DMFT) has been proved to be a very useful and effective tooP^"^^, which has 
made significantly progress as in the study of metal-insulator phase transition. The DMFT is exact when 
investigating the strongly correlated system in the infinite- dimension, in which the self- energy is independent 
of momentum. However, in low- dimensional systems, the quantum fluctuation and short range correlations play 
important roles, which are ignored in DMFT. The cellular DMFT (CDMFT), as a cluster extension of DMFT, 
effectually incorporates the spatial correlations by mapping the many-body problem into finite clusters embed- 
ding in a self-consistent media^^"^^. In two-dimensional systems CDMFT is more precise than DMFT, and is more 
effective to investigate the phase transition in low and multi-component systems^ 

In this report, we investigate the finite temperature metal-insulator and magnetic phase transition in strongly 
correlated bilayer honeycomb lattice (BHL) by combining the cellular dynamical mean-field theory (CDMFT) 
with continue-time quantum Monte Carlo (CTQMC) method^^. A phase transition from paramagnetic phase to 
anti-ferromagnetic phase occurs by investigating the magnetization. In a proper value of inter-layer hopping, a 
novel layer anti-ferromagnetic phase emerges. The layer anti-ferromagnetic phase can be broken by the increas- 
ing inter-layer hopping, and a paramagnetic phase can be found. Moreover, the nonlocal inter-layer hopping 
plays an important role on localizing the free election and modif)^ing the spatial distribution of the electron in 
lattice sites, especially in dimer sites. For instance, the dimer sites of our system are double occupied, and the non- 
dimer sites are single occupied. We have also presented the DOS, double occupancy, and fermi surface below, 
which can be directly detected in future experiments. 

Results 

The strongly correlated bilayer honeycomb lattice. Fig. 1(a) shows the lattice structure of BHL. Ai and Bi 
denotes the sites on top-layer, while A2, B2 signs the bottom-layer sites. The inter-layer band connects Ai (top) 
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Figure 1 | The structure of bilayer honeycomb lattice and its qualities in 
the non- interacting limit, (a): Bernal stacking of the bilayer honeycomb in 
real-space with intra- and inter-layer hopping t and ti between the 
sublattice Ai, Bi on top-layer and A2, B2 on bottom-layer. The black arrows 
Ui and a2 are the lattice vectors, (b): The cellular of our bilayer system in 
cellular dynamical mean-field theory (CDMFT). The red solid line 
containing sites 0, 1, 2, 3, 4, 5 belongs to top-layer and blue dotted line with 
sites 6, 7, 8, 9, 10, 11 belongs to bottom-layer, (c): Reciprocal lattice of 
bilayer honeycomb lattice. The thick red line shows the first Brillouin zone. 
The r, Ky M and K' points denote high symmetry points in first Brillouin 
zone, (d): Density of states of our system for A1/A2 and B1/B2 sites when U 
= 0. 

and A2 (bottom) sites. The first Brillouin zone can be found in 
Fig. 1(c), in which the F, iC, M and K' shows the points with 
different asymmetry in the A;-space. The non-interaction density of 
states (DOS) for different sites when ti = t is shown in Fig. 1(d), 
which is different from the mono-layer honeycomb lattice^^. The 
low- energy dispersion is quadratic, which is linear in the mono- 
layer case. There is no band gap between the conduction and 
valence bands in BHL. In this report, we investigate the phase 
transition on the half-filling BHL. We consider the standard 
Hubbard model: 



H = 



'i(joc^i(j{l-a) ■ 



where t is the intra-layer nearest-neighbor hopping energy, U is the 
Coulomb interaction, ti denotes the inter-layer hopping energy, fi is 
the chemical potential, which is adjusted to keep the system at half- 
filling, c]^^ and Cifjoc denote the creation and annihilation operators of 
fermionic atoms on site / with spin cr respectively, and a shows the 
layer parameter (a = 0 means the top -layer, and a = 1 denotes the 
bottom-layer), niaa = cl^oc^iaa corresponds to the density operator. 
Here we set ^ = 1, which is also used as the energy unit in this report. 

The metal-insulator phase transition. The ti — (7 phase diagram is 
shown in Fig. 2. The phase boundary of the metal and insulator is 
shown by the red solid line with square plots. When = 0 the critical 
point of metal-insulator phase transition is 11^ = 4.7, meaning that 
BHL degenerates into a mono-layer one. This result is the same as 
Ref 16. In order to study the magnetic phase transition, we define the 
staggered magnetization as M = S/S^n(/)(< > — < >), where 
/ denotes the lattice index, sgn{i) = 1 for A 1 and A2 sites and sgn{i) = 




Figure 2 | Metal- insulator phase diagram in fixed Tor ti. Phase diagram 
as a function of inter-layer hopping ti and interaction (7 at T = 0.1. The red 
solid line shows the boundary of metal-insulator phase transition, and the 
blue solid line denotes the magnetic phase transition, which divides the 
phase into paramagnetic metal (PM: M = 0 and AE =0), paramagnetic 
insulator (PI: M = 0 and AE 7^ 0), anti-ferromagnetic metal (AFM: Mt^O 
and AE = 0) and anti-ferromagnetic insulator (AFI: M 7^ 0 and AE 7^ 0). 
Inset: Phase diagram as a function of temperature T and interaction U at 
fixed ti. 

— 1 for Bi and B2 sites. The phase boundary of non-magnetic state 
and magnetic state is shown by blue solid lines with circle points. As 
shown in Fig. 2, when 2 < ti < 4.8, anti-ferromagnet can be found. 
When U is increased, a phase transition from gapless anti- 
ferromagnetic metal (AFM), where AE = 0, M 7^ 0 and AE here is 
single particle gap, to gapped anti-ferromagnetic insulator (AFI), 
where AE 7^ 0, M 7^ 0, occurs. When < 2, a paramagnetic metal 
(PM), where AE = 0, M = 0, is formed, and a PM-AFM-AFI 
transition happens when U is increased. In the value of 0 < < 
0.2, when U > 4.7, a small region, which is named as paramagnetic 
insulator (PI) is found with AE 7^ 0 and M = 0. 

Double occupancy is always used to measure the localization of the 
electrons directly, which is an important parameter to observe the 
metal-insulator phase transition^^. In this report, we investigate the 

double occupancy Docc = ^F / dU = - HfKni^nii > as a function of 




Figure 3 | The evolution of double occupancy Docc The double 
occupancy as a function of inter-layer hopping ti for different interaction 
17 at temperature T = 0.1. The dimer sites tend to be double occupied 
however non- dimer sites tend to be single occupied, with increasing ti. 
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Figure 4 | The density of states, (a) - (c), (e) - (g): The density of states as a function of frequency co for different inter-layer hopping ti at temperature T 
= 0.1. (a), (e): The metallic phase at U = 2.5 for A sites and B sites, (b): For A sites, as we increase ti, system undergoes a phase transition from 
metal (at U = 3.5 and ti = 1.8) to insulator (at ti = 3.2). Single particle excitation gap opens at about ti = 3.2. (f): For B sites, same tendency occurs and 
the critical point is also at about ti = 3.2. (c), (g): The insulating phase at 17 = 6.0. There is a visible single particle excitation gap around the Fermi energy, 
(d), (h): The density of states as a function of co for different 17 with fixed ti. 



inter-layer hopping ti for various interaction U [see in Fig. 3], where 
F is the free energy. We use D^^^ to describe the D^cc for Ai and A2 
sites (hollow circles in Fig. 3), and Dbb the D^cc for ^1 and B2 sites 
(solid circles in Fig. 3). When ti = 0, D^cc decrease when U is 
increased showing the localization is enforced by the interaction. 
When ti 7^ t, D^^ is different with Dbb- It is found that when ti is 
increased, D^^ increases while Dbb decreases. This suggests the itin- 
erancy of electrons in dimer sites is enhanced due to the increasing 
inter-layer hopping. However, be different with D^^, Dbb decreases 
while ti is increased. These results suggest the intra-hopping between 
dimer and non- dimer sites is weaken due to the forming of spin- 
polarized electrons in dimer sites. 

We derive DOS from the imaginary time Green's function G(t) to 
observe the single particle spectral for different interaction U and ti 
by the maximum entropy method^^'^^. Fig. 4(a) and Fig. 4(e) shows 
the DOS for different inter-layer hopping when U = 2.5, T = 0.1 for 
A sites and B sites. It is found that, in weak interaction, the inter-layer 
hopping ti does not affect the metallic properties of BHL. In Fig. 4(b) 
and Fig. 4(f), we can find that the system keeps at a metallic state 
when ti = 1.8, and an obvious pseudo-gap is found when ti = 2.2. A 
metal-insulator transition happens when ti = 3.2 for U = 3.5 and T 
= 0.1 for both A and B sites. When interaction being large, such as U 
= 6.0 [see Fig. 4(c) and Fig. 4(g)], the system stays at an insulating 
state, which is insensitive of ti. Fig. 4(d) and Fig. 4(h) show the DOS 
for different (7 when = 1.0 and T = 0.1. An opened gap is found 
when the interaction is increased, which indicates a phase transition 
from metal to insulator. These results mean that the inter-layer hop- 
ping and Coulomb interaction play same roles for the metal-insu- 
lator phase transition. 

In this report, the Fig. 5 shows the spectral function A(k, co = 0) 
near the Fermi surface for various U and ti. A(k, oj) is defined as: 

A(A;,(X))~- -lim ImGk(o) + iO) 



Gk{ico) 



1 



/(X> — e^ — S^(/co) 



where ico is the Fermionic Matsubara frequency, denotes the dis- 
persion relation, T^ki^co) corresponds to the /c-dependent self-energy, 
k is the wave vector in the first Brillouin zone. Figs. 5 (a) -(c) show the 



A{ky CO = 0) for different ti when U = 3.5. When ti < 3.2, we can 
find that the peaks of A(k, oj = 0) distribute in the K and K' 
points [see Fig. 1(c)] in /c-space, displaying a metallic behavior 
[see Figs. 5(a), (b)]. When ti = 4.0, the A{k, cd = 0) ^ 0 [see 
Fig. 5(c)], meaning that a gap can be found near the Fermi energy. 
This gapped behavior shows the system is an insulator. Similar 
behavior is found when ti = 1.2 and U is increased [see Figs. 5(d) 
- (f)]. It means both spin -polarization and strongly correlated 
interaction can induce the localization of electrons, driving the 
system from metal to insulator. 

The magnetic phase transition and a novel layer anti-ferromagnetic 
phase. In strongly correlated BHL, charge imbalance between the 
two sublattices sites causes different kinds of magnetic spontaneous 
symmetry breaking, dividing the sites into dimer sites and non- 
dimer sites. In order to study the magnetic order, we use a magne- 
tic order parameter defined as = ^ ( < n/| > — < > ) , 

where a denotes the sublattice Aj, A2, Bi and B2, is the lattice index 
for sublattice a, and N^, means the site number for sublattice a [see in 
Fig. 1(a) and Fig. 1(b)]. Parameter < > indicates the electron 
density in lattice site / with spin index a (we set magnetization of Aj 
as positive sign). Fig. 6 shows the evolution of nf as a function of ti 
for U = 3.8. It is found that when U = 3.8, m is zero for ti < 1.0, 
denoting that no magnetic order is formed for weak ti. A magnetic 
state with anti-ferromagnetic order is found when ti > 1.0. When ti 
> 4.4, the m^^ and m^^ decrease to zero, while m^^ are still 
finite, meaning a phase with non-magnetic order and magnetic order 
coexisting. This novel phase is called layer anti-ferromagnetic state, 
in which the dimer sites are non-magnetic and non- dimer sites are 
magnetic. The sketches of the possible magnetic order existing in 
BHL is shown in Figs. 7(a) - (c). 

Finally, the phase diagram of magnetization about ti and L^for T = 
0.1 is shown in Fig. 7. A phase transition from paramagnet to anti- 
ferromagnet is found when U is increased. A layer anti-ferromag- 
netic phase is found for large ti, in which the magnetic of dimer sites 
keeps zero while non-dimer sites is nonzero. When U > 3.7, a phase 
transition from layer anti- ferromagnetic state to paramagnetic state 
is found. 
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Figure 5 | The evolution of spectral function near Fermi surface. The spectral function A{k, w = 0) near the Fermi surface for different inter-layer 
hopping at r= 0.1 and U= 3.5: (a) h = 1.0, (b) ti = 1.8, (c) fi = 4.0. The A(fc, co = 0)whenr= 0.1 and /"i = 1.2 for different interaction 17: (d) U= 2.0, 
(e) U= 2.8, (f) 17= 4.5. 



Experimental protocol. We propose an experiment setup to 
investigate the phase transition in strongly correlated bilayer 
honeycomb lattice (BHL). The atoms can be produced as a 
pure fermion condensate by evaporative cooling^^, which provides 
two hyperfme states \F,mp) = |9/2, -9/2) = | j) and \F,mp) = |9/ 
2, — 7/2) = I ]y^. Three standing- wave laser beams are used to form 
the honeycomb lattice, and two extra laser beams along the z 
direction suppress the tunneling between layers^^. The potential of 
optical lattice is given by Vh(x, y) = VoS^=i,23 sin^[k(x cos Oj+y sin 
0j) + 7i/2], where = 7r/3, O2 = 27i/3, 9^ = 0. Then, we use another 
three standing-wave laser beams with a 27i/3 angle between each 
other to form triangular lattice. The potential is given by 

Vt{x,y) = Vq 1^3 +4 cos{kxx/2) cos (y3kyy /2) + 2 cos (a/3A;^)/) ] . K 

and ky are the two components of the wave vector k = 2nlX in these 
two types of lattices, where X = 73S nm is wavelength of the laser, 
and Vq is given in recoil energy £^ = h^k^llm. Inserting the triangular 
lattice between two layers of honeycomb lattice, the Bernal stacking 
BHL with trapped^^K atoms can be formed^^'^^. In BHL the intra-layer 

hopping t= {4./ y/^)Elf'^vl'^exp^-2(VQ/Erf'^'^ is adjusted by the 

periodic potential of laser beam and ti can be tuned by changing the 
wavelength of laser beam in z direction. The on-site interaction 

U= \/S/ nkasEr{Vo / ErY^"^ is determined by the s-wave scattering 



length which can be tuned by Feshbach resonance. The 
temperature can be extracted from the time-of-flight images^ \ 

It should be mentioned that, we can detect the numbers of double 
occupied sites to confirm whether the metal-Mott insulator trans- 
ition happens. Firstly, we have to increase the depth of the optical 
lattice to prevent further tunneling of atoms. Next, the energy of the 
atoms on doubly occupied sites is shifted by approaching a Feshbach 
resonance. Then the one spin component of atoms on double occu- 
pied sites is transferred to a new magnetic sublevel by radio - 
frequency pulse method. Finally, the double occupancy can be 
deduced by the absorption images^^'^^ 

To get Fermi surface in experimental, we ramp down the optical 
lattice slowly enough first, and the atoms stay adiabatically in the 
lowest band while quasi-momentum is approximately conserved. 
Then the lattice potential is lowered to zero rapidly, by switching 
off the confining potential and the atoms can expand ballistically for 
several milliseconds. The Fermi surface can be obtained by a absorp- 
tion image^^'^^. 

Discussion 

In this work, we have investigated the metal-insulator transition and 
magnetic phase transition in strongly correlated bilayer honeycomb 
lattice using cellular dynamical mean-field theory (CDMFT) 
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Figure 6 | The evolution of the magnetic order parameter m. The 

evolution of magnetic order parameter m at T = 0.1 and U = 3.8. A 
paramagnetic phase is found in weak ti. For ti > 1.0, the magnetic 
parameter m is nonzero and has opposite sign between A1M2 sites and BJ 
B2 sites. The system goes into anti-ferromagnetic phase. At large ti the 
magnetization of A1/A2 sites are more easily decreasing to zero while B1/B2 
sites are still nonzero. The system is layer anti-ferromagnetic phase. Single 
particle excitation gap AE denoted by the dark green solid line and orange 
solid line, divide the phase into paramagnetic metal (PM), anti- 
ferromagnetic metal (AFM), anti-ferromagnetic insulator (AFI) and layer 
anti-ferromagnetic insulator (LAFI). 

combining with continue-time quantum Monte Carlo (CTQMC) 
method. In lower temperature case we map the phase diagram as a 
function of interaction U, inter-layer hopping ti and magnetization 
m. It shows that the inter-layer hopping affects the electrons to form 
spin-polarized electrons, and an insulating state is induced. A layer 
anti-ferromagnetic phase is found at large ti, in which the magnet- 
ization of dimer sites is zero while non-dimer keeps finite value. 
Therefore, the inter-layer hopping ti plays an important role to form 
a singular magnetic spontaneous symmetry breaking phase. Our 
study may provide a helpful step for understanding the interaction 
and inter-layer hopping driven metal-insulator transition, the exotic 
magnetic order with asymmetry breaking and the possible magnetic- 
nonmagnetic order coexisting state. 

Methods 

The cellular dynamical mean-field theory. We combine the cellular dynamical 
mean-field theory (CDMFT) with continuous time quantum Monte Carlo (CTQMC) 
method to determine the metal-insulator transition and magnetic phase transition in 
the strongly correlated bilayer honeycomb lattice. In low- dimensional systems, 
quantum fluctuations are much stronger than the higher dimensions. The nonlocal 
effect is much important in this case. Dynamical mean-field theory ignoring the 
nonlocal correlations leads lots of errors in calculation. Therefore, we use CDMFT, as 
the advanced method in our work. We map the original lattice onto a 12-site effective 
cluster embedded in a self-consistent bath field [see Fig. 1 (b)] . Starting with a guessing 
self-energy (which is independent of momentum^^), we can get the Weiss field 
Goiioj) obtained by the coarse-grained Dyson equation: 




Layer Anti-Ferromagnet 
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Anti-Ferromagnet 
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U 

Figure 7 | The phase diagram of magnetic phase transition. In weak 
interaction [/and weak inter-layer hopping ti, the system is paramagnet. 
When we increase 17, the system undergoes a magnetic phase transition to 
anti-ferromagnet. When we increase ti, the magnetization of A1/A2 sites 
decrease to zero while B1/B2 sites stay finite. The system becomes layer anti- 
ferromagnet. In the region where ti > 5.0, the system returns to 
paramagnetic phase. 



i(o — t{k) — lL{i(o) 



(1) 



where ico is Fermionic Matsubara frequency, /i is the chemical potential, k is in the 
reduced Brillouin zone of the super-lattice, and t{k) is hopping matrix for the super- 
lattice. The form of t{k) is: 
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and ai, a2 are real lattice vectors as 



} (2) 



where ( 

in Fig. 1(a). The cluster Green's function G{ico) can be gotten by the impurity solver. 
In our work, we use the numerically exact CTQMC simulation as impurity solver and 
take 5 X 10*^ QMC sweeps for each CDMFT loop'*^. The new self-energ)^ Z(/a») is 
recalculated by the Dyson equation: 

i:{ico) = Q- ^ (ico) - G" ^ (ico). (3) 

This iterative loop repeated until self-energy is converged. 

The CTQMC method as impurity solver can be taken as follows. We start the 
procedure at partition function, which can be written as: 



Z = Tre 



{T)dT 



(4) 
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where T^. is time-ordering operator, Hi (t) = e^^^Hie is Hi in the interaction 
picture, and Zq = Tre~^^° is a partition function for the unperturbed term. Putting Hi 
= UT^iHiirii-^ in Eq. 4, the partition function is 

^ = ^0 E^^^j- • • l^^i • • • M^Mri) ■ ■ ■ n^{n))^{T,n^{r,) ■ ■ ■ n^{r,)) ^. (5) 

Here (>o indicates a theromdynamic average with respect to e~^"\ Using Wick's 
theorem, for each order in k, {TT^Uairi) ■ ■ ■ na(rk))Q((^ = ],i) can be written as 
determinant detD{k): 



m- 



fG«(ri,ri) G«(ri,r2) 
G\r2,r,) G'{r2,r2) 



G\r2jk) 



(6) 



where G° is non-interacting Green's function. There is no spin index in D{k) for the 
determinants of spin-un and -down being equivalent. Like classical Monte Carlo, by 
integrand of Eq. 5, we can get the weight of order k 

Wk = {- dxUfdetD^ {k)detDi (k), (7) 

where St = jS/L is slice of imaginary time. We can get the standard Metropolis 
acceptance ratio R of adding vertex by the detailed balance condition: 

j^WkPk^k+i = j^Wk+iPk+i^k, (8) 



Pk+i^k k-\-l\ detD^{k)detDi{k) 



(9) 



Here P^-^k+i is the probability to increase the order from A; to A; + 1 {Pk+i->k the 
probability to decrease the order from k + 1 to /c), is probability to choose a 

position in time and space for vertex you intend to add while ^^py probability to 

choose one vertex you intend to remove of from the existing k + I noes. To calculate 
the ratio R, we have to deal with the function detD{k + l)/detD{k). 

detD{k+ 1) /detD{k) = detiI^{D{k^ 1) -D{k))M{k))=X, (10) 
M(;(k) = D~ ^ (k), we can easily get the value of /I in matrix form: 



det 



0 0 

\G\n+uri)M{k)., G\n+uri)M{k)^, 



1 G«(rfc,rfc+i) 
G\rk+uri)M{k)., G\n+urk+i) J 



(11) 



= G\rk+urK+i)-G\rk+uri)M{k).^G%rj,n+,)=L 

Then it is easy to obtain the update M for the order /c + 1 by numerical method: 



M(A:+1) 



/ 



\-r'Rk+i 



M'i,j 



-^l,k+\X 

-Lk,k+iX~ 



-r 



(12) 



where the factor of the matrix is M'i^j = M{k)ij + Li^k+iX ^Rk+i,j, Ri,j = G%i, l)M{k)ij 
and Li^j = M{k)ijG"{l, j). For the step k — 1, we can also get the radio R and update 
formulas of M{k - 1): 



^_ k fdetD^(k-l)detD^(k-l) 
^~ UPN \ detD^{k)detDi{k) 



Mi^j(k - 1 ) - Mi^j (k) - Mij {k)Mi^j (k) /Mij (k) . 



(13) 
(14) 



Using the update formula for M, the Green's function can be obtained both in 
imaginary time and at Matsubara frequencies: 



G(t - t') = G° (t - t') - G° (t - Ti)MijG^ {zj - t') , 



G{icQ) = Go{icQ)-Go{icQ) 
Here Go(/a;) is a bare Green's function. 



Go ((to). 



(15) 
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